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We propose a new implementation of the replica-exchange method (REM) in which replicas follow 
a pre-planned route in temperature space instead of a random walk. Our method satisfies the 
detailed balance condition in the proposed route. The method forces tunneling events between 
the highest and lowest temperatures to happen with an almost constant period. The number of 
tunneling counts is proportional to that of the random-walk REM multiplied by the square root of 
moving distance in temperature space. We applied this new implementation to two kinds of REM 
and compared the results with those of the conventional random-walk REM. The test system was 
a two-dimensional Ising model, and our new method reproduced the results of the conventional 
random-walk REM and could adjust the tunneling counts by two times or more than that of the 
random-walk REM by replica-exchange attempt frequency. 
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INTRODUCTION 


The enhancement of configurational sampling ensures accuracy and efficiency of molecular dynamics (MD) or Monte 
Carlo (MC) simulations. Replica-exchange method (REM) [1-4] (or parallel tempering) is a popular way to improve 
efficiency of configurational sampling (for reviews, see, e.g., Refs.[5, 6]). This method achieves a random walk in 
temperature space and allows the system to overcome energy barriers between local-minimum free energy states. The 
number of tunneling events [7, 8] is defined to be the number of times where the simulation visits from the lowest 
temperature through the highest temperature and back to the lowest temperature (hereafter, we refer to this number 
as tunneling counts). The tunneling count is a measure of configurational sampling efficiency: the more tunneling 
counts are obtained with a fixed number of simulation steps, the better the sampling is. However, as the number of 
replicas increases, tunneling events require more steps by random walk. This is a difficulty for large systems because 
they require many replicas in general. 

In the Metropolis criterion with pseudo random numbers, randomness prevents the design of a trajectory in tem¬ 
perature space. Moreover, the moving distance by a random walk is proportional to the square root of the number 
of trials. Therefore, previous improvements of tunneling count and efficiency in the REM was fulfilled by approaches 
such as the temperature selection[9-12], attempts of non-neighboring pairs, and increase of acceptance probabilities. 
An example is all-paring method, which tries all combinations of temperature pairs[13, 14]. Another example is Gibbs 
sampling heat-bath replica-exchange method [15]. Other examples are to employ the global balance condition such 
as replica permutation method[16] with Suwa-Todo algorithm[17] and all possible pair exchange[18]. These methods 
reported higher tunneling counts than the conventional REM. However, a production of an equilibrium state without 
random walk leads to another possibility of an improved REM, which enables a simulation to determine a trajectory 
of replicas in a temperature space. In fact, some of special manipulation of state changes without the detailed balance 
condition shows a rapid convergence for thermal equilibrium[19, 20]. Based on Chaotic Boltzmann machine[21, 22], 
we have also proposed the deterministic replica-exchange method (DETREM)[23], which performs replica exchange 
based on a differential equation without pseudo random numbers. This method produced the same efficiency as the 
conventional REM. 

Recently, Spill et al. showed that a planned route trip in only one randomly selected replica during a simulation 
showed an improvement compared with the conventional REM [24]. We can generalize this idea so that planned 
route for all replicas will realize the maximum efficiency with wide configurational sampling, and we here propose the 
designed-walk replica-exchange method (DEWREM) by even-odd sequential replica exchange. 
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METHODS 


We now give the details of our methods. We prepare M non-interacting replicas at M different temperatures. Let 
the label i (=1, • • •, M) stand for the replica index and label m (=1, • • •, M) for the temperature index. We represent 
the state of the entire system of M replicas by X = w here itl = {qM,p[*)} m are the set of 

coordinates and momenta pW of particles in replica i (at temperature T m ). The probability weight factor for state 
X is given by a product of Boltzmann factors: 

M 

W REM (X) = JJexp[ -/3 m (i)H(qW,pW)\, (1) 

i =1 


where /J m (= l/k^,T m ) is the inverse temperature and H(q,p) is the Hamiltonian of the system. We consider exchanging 
a pair of replicas i and j corresponding to temperatures T m and T„, respectively: 
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wherein = |gW,pH'| , and pW = [4]- 

Here, the transition probability ui(X —> X') of Metropolis criterion for replica exchange is given by 


where 


uj{X —> X') = min 


/ Wrkm{X') \ 
V ’ Wrem(-X’) / 


= min(l, exp(—A)), 


(3) 


A = A m , n = (/?„ - p m ){E{qW) - E{q®)). 


(4) 


Because each replica visits various temperatures followed by the transition probability of Metropolis algorithm, REM 
performs a random walk in temperature space. 

We now review two REMs, which are based on random walks in temperature space. Without loss of generality, we 
can assume that M is an even integer and that Tf < T 2 < ■ ■ ■ < Tm- The conventional REM[1-4] is performed by 
repeating the following two steps: 


1. We perform a conventional MD or MC simulation of replica * (= 1, • ■ ■ , M) at temperature T m ( m = 1, ■ ■ ■ , M) 
simultaneously and independently for short steps. 

2. Pairs of exchange attempts are selected in replica pairs with neighboring temperatures, for example, for the odd 
pairs (Ti,T 2 ), (T 3 ,T 4 ),---, (T M -i,T m ) or even pairs (T 2 ,T 3 ), (T 4 , T 5 ),-• •, (T M _ 2 ,T M _ 1 ). 


All the replica pairs thus selected are attempted to be exchanged according to the Metropolis transition probability 
in Eqs. (3) and (4) with n = m + 1. 

We repeat Steps 1 and 2 until the end of the simulation. The canonical ensemble at any temperature is reconstructed 
by reweighting techniques [25-27]. 

We next present the deterministic replica-exchange method (DETREM)[23]. Only Step 2 is different from the 
conventional REM. At first, we introduce an internal state y m ,n as an index of a pair of replicas i and j at temperatures 
T m and T n , and consider the following differential equation: 

dym,m +1 1 / c \ 

"TT — CT m - ■ 77 r, (5] 

at 1 + exp(A mjm+ i) 

where t is a virtual time, A mim+ i is the same as in Eq. (4) with n = m + 1, and the signature a m of the pair of 
(T m , T m+ i) changes to 1 or —1 to control the signature of the change of y m which monotonically increases or decreases. 
In Step 2, instead of applying the Metropolis criterion in Eqs. (3) and (4), we solve the differential equation in Eq. 
(5) for the internal states y m ,m+i G {—1,1} for (T m ,T m+ i), where the total number of internal states is M- 1 with the 
following pairs: (1,2), (2,3), • ■ •, ( M-1,M) for the random-walk DETREM and the pairs: (1,2), (3,4), • • •, (M-1,M) 
and (2,3), (4,5), • • •, (M-2,M-1) for designed-walk REM. The replica exchange is done as follows[23]: 

If updated ym,m -\-1 ^ dil, then [T rni T rn j r \ s j y (T m _|_ 4 ,T m ), 

ym,m-\- 1 ^ 1 “F I 5 <-Tl- 
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For the random-walk DETREM, if y m ,m +1 performs exchanges, y m +i,m +2 is not time evolved and y m + 2 ,m +3 is evolved 
to avoid the leap exchange of temperature such as from T m to T m+2 . 

Finally, the designed temperature walk can be implemented to both conventional REM and DETREM (and other 
REMs) as follows. Namely, the designed-walk replica-exchange method (DEWREM) is performed by repeating the 
following steps. 

1. We perform a conventional MD or MC simulation of replica i (= 1, ■ • • , M) at temperature T m (m = 1, ■ ■ • , M) 
simultaneously and independently for short steps. 

2. Replica exchange is attempted for all the odd pairs (Ti,T 2 ), ( 13 , 14 ), • • •, (7m-i,1m). 

3. Repeat Steps 1 and 2 until all odd pairs perform replica exchange exactly once. Namely, once a pair is exchanged, 
the exchanged pair stops exchange attempts and keep performing the simulation in Step 1 with the new temperatures. 
Replica exchange attempt in Step 2 is repeated until all the other odd pairs finish exchanges. 

4-6. Repeat Steps 1-3 where the odd pairs in Steps 2 and 3 are now replaced by the even pairs (T 2 ,T 3 ), ( 74 , 15 ), 

■ ■ •, (7m-2,7m_i). 

7. The cycle of Steps 1 to 6 is repeated until the number of cycles is M, which is equal to the tunneling count and 
all replicas have the initial temperatures. 

8 . Begin the above cycle of Steps 1-7 with Steps 1 to 3 and Steps 4 to 6 interchanged. 

These eight steps are repeated until the end of the simulation. 

The schematic picture of this procedure is shown in Fig. 1. We remark that Step 8 , namely, reversing the cycle of 
Steps 1-3 and 4-6, is necessary for the detailed balance condition, because the entering states are the same as leaving 
states. For example, the state ( 2 ^, 2 !^, 2 !^, x 1 ^], 2 ^, xj?) is reached from only two states (x^, x^, xj^, x^, 2 ^, 2 ^), 
( 2 ^, xjj 2 ^, 2 ^, 2 g 4 ^, 2 ^, Xg 6 ^) only and makes transition to the two states as shown in Fig. 1, where Xm is the state 
of replica i at temperature T m . This exchange procedure satisfies the detailed balance condition for replica and 
temperature pair because the trials of exchange pair 


7 ( i(m) —> i(m + l))w ((x^ m) , x^ 1] ) -+ ( 2^1 , x^ ro+1) )) 

= 7(i(m + l) -+ i(m))w ((^+i,a;^ m+1) ) -+ (x^ m) ,x^+ 1) )) 


( 6 ) 


is equal in the route as is shown in Fig. 1, where 7 ( 1 ( 777 ) —► 1 ( 777 + 1 )) is the selected probability of the exchange attempt. 
Note that for exactly the same conditional probability of odd pair replica exchange, (Ti,T 2 ), ( 73 , 74 ), • • •, (Tm-i,Tm), 
A is J2k=i (P 2 k-i — p 2 k)(E 2 k — Ej 2 k—i)i where k = 1, • • ■ , M/2. However, because each replica is non-interacting with 
other replicas, waiting for other exchanges does not influence the transitions of others. 

This sequential exchange achieves one tunneling count when M cycles for each replica are finished. In theory, the 
estimated ratio of tunneling count between the odd-even sequential exchange and the conventional random walk is 
given by 
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where Aerial is the number of exchange attempts, -P c ^5^ tion is the correction for waiting for all the replica exchanges 
in Steps 3 and 6, and -P c ^^ ectiorl is the correction for the deviation of random-walk probability from the value 1/2. 


RESULTS 

In order to test the effectiveness of the present methods, we applied them to the 2-dimensional Ising model. The 
lattice size in a square lattice was 128 (hence, the number of spins was N = 128 2 = 16384). We have performed 
conventional random-walk simulation and designed-walk (DEWREM) simulation of both conventional REM and 
DETREM. We have also performed a mixed random-walk and designed-walk simulation of DETREM, where we 
repeated the two walks alternately. The total number of replicas M was 40 and the temperatures were 1.50, 1.55, 
1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.94, 1.98, 2.01, 2.04, 2.07, 2.10, 2.13, 2.16, 2.19, 2.22, 2.25, 2.28, 2.31, 2.34, 
2.358, 2.368, 2.38, 2.40, 2.42, 2.44, 2.47, 2.51, 2.57, 2.63 ,2.69, 2.75, 2.82 ,2.90, 3.00, 3.10, and 3.15. Boltzmann 
constant fcs and coupling constant J were set to 1. Thus, /3 = l/A+T 1 = 1/T = /3*, and the potential energy is given 
by E( s) = — j> s i s ji where Sj = ±1, and the summation is taken over all the nearest-neighbor pairs in the square 
lattice. 
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For the conventional random-walk REM and DETREM, replica-exchange attempt was made every 1 MC step. One 
MC step here consists of one Metropolis update of spins. The total number of MC steps for all the simulations was 
100,000,000. To integrate Eq. (5), we used the fourth-order Runge-Kutta method with virtual time step dt = 1. 
For DEWREM, replica-exchange attempt was made every 10, 20, 50, 100, and 150 MC steps in the conventional 
REM simulations and every 10, 20, 50, and 100 MC steps in the DETREM simulations (see Table I). The mixed-walk 
simulation was performed in which after 4 M{= 160) even-odd or odd-even cycles of designed-walk simulations (replica- 
exchange attempt was made at every 20 MC steps) were performed, 200,000 MC steps (which roughly corresponds 
to 2 M cycles) of random-walk simulations (replica-exchange attempt was made at every 1 MC step) were performed, 
and then this procedure was repeated. For reweighting analyses[25-28], the total of 10,000 spin state data were taken 
with a fixed interval of 1,000 MC steps at each temperature from the REM simulations. 

Table I lists the mean tunneling counts in temperature space and energy space per replica for each method, which 
is the number of times where the replicas visit from the lowest temperature through the highest temperature and back 
to the lowest, and the mean energy at the lowest temperature through the mean energy at the highest temperature 
and back to the mean energy at the lowest temperature during the simulation, respectively. The mean tunneling 
counts in temperature space per replica of the designed-walk simulations at every 100 MC attempts were about four 
times larger than random-walk Metropolis REM and twice larger than random-walk DETEM. On the other hand, 
the mean tunneling counts in energy space per replica of the designed-walk simulations at every 100 MC attempts 
were about six times larger than random-walk Metropolis REM and twice larger than random-walk DETEM. These 
large tunneling counts imply that in designed-walk method all replicas traversed more efficiently in temperature 
space, and our design to maximize the tunneling counts for all replicas without random walks was successful. For the 
mixed-walk simulation, the maximum tunneling count was about twice larger than that of random-walk DETREM. 
The mean tunneling count was almost the same as that of designed-walk DETREM. We next examine physical 
quantities obtained from the designed-walk simulations with various replica-exchange attempt frequencies and mixed 
walk simulations and compare them to those from the conventional random-walk simulations. Fig. 2(a) and Fig. 2(b) 
compare the canonical probability distributions of energy density e = E/N at four temperatures as a function of T and 
the average total energy density e as a function of T obtained from the random-walk and designed-walk simulations 
of REM and DETREM, for DETREM including the mixed-walk simulation. Most of the probability distributions 
of energy density are the same. However, in the DEWREM simulations, the probability distributions with high 
frequency replica-exchange attempts at T — 2.25 near the exact critical temperature T c = 2.269 are deformed slightly 
compared to the results of the random-walk simulation. This can also confirm the results of heat capacity as shown 
in Fig. 3. Fig. 2(c) and Fig. 2(d) show the average total energy density e = E/N as a function of T from the same 
simulations. They were obtained by the reweighting techniques[25-28]. The average total energy densities in all the 
simulations are the same. 

Fig. 3(a) and Fig. 3(b) show the specific heat C as a function of T during the conventional REM simulations and the 
DETREM simulations, respectively. This shows that designed-simulation with shorter replica-exchange interval such 
as every 10 and 20 MC steps underestimated the heat capacity near the critical temperature although the transition 
point is sufficiently similar to the exact critical temperature at T c = 2.269. As the intervals of replica-exchange 
attempts are longer, the accuracy of heat capacity is higher. Moreover, the combination of the random-walk and 
designed walk also increased the accuracy. This suggests that the designed-walk replica-exchange attempts caused 
correlation between replicas. The correlation seems to be very strong near the critical temperature. As a result, 
the heat capacity is underestimated slightly. However, these physical quantities show that mixed-walk simulation can 
increase the accuracy of results and the number of tunneling counts and DEWREM simulation is suited for simulations 
with longer replica-exchange time intervals between replica-exchange attempts. 


CONCLUSIONS 

In this article, we proposed an algorithm by a designed temperature route for replica-exchange methods. The aim for 
developing this method was to maximize the tunneling counts. We reproduced the results of conventional REM. Our 
designed-walk simulations showed that the tunneling counts increased by at least twice more than the random-walk 
simulations. This method also showed that in REM replicas have no interaction among replicas but replica-exchange 
trials with artificial orders cause a correlation between replicas. By decreasing the correlations introducing mixed-walk 
simulation, using longer- interval replica-exchange trials, or other paring of replicas, this method will be more efficient. 

In a future work, we will give the formulation for designed routes for multidimensional replica-exchange method[29], 
which is a power tool in massive parallel computing [30-32]. This will be useful for creating another design route for 
non-neighboring update for exchanges. 
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FIG. 1. An schematic picture of time series of temperature indices in DEWREM with 6 replicas. The left cycle begins with the 
temperature exchange of odd index pairs (Ti,T 2 ), ( Tz,T ±), and (T^,Tq), then tries with even pairs {T 2 ,Tz) and (Ti^Ts). The 
right cycle begins with even pairs and next tries odd pairs. They are the reverse cycles of each other and their combination 
satisfies the detailed balance condition of replica exchange. 


TABLE I. The mean number of tunneling counts per replica. 


TC 

Random walk 

Designed walk 

Mixed 

T TC 

Met 

DETREM 

Met 

DETREM 

DETREM 

Interval 

1 

100 

1 

100 

20 

50 

100 

150 

20 

50 

100 

1 & 20 

Mean 

173 

37 

178 

58 

292 

197 

131 

99 

231 

144 

93 

293 

± SD 

10 

3 

9 

5 

56 

41 

27 

21 

48 

30 

20 

6 

E TC 

Met 

DETREM 

Met 

DETREM 

DETREM 

Interval 

1 

100 

1 

100 

20 

50 

100 

150 

20 

50 

100 

1 & 20 

Mean 

55 

13 

55 

34 

80 

79 

77 

70 

81 

79 

69 

75 

± SD 

5 

2 

5 

4 

8 

6 

7 

5 

6 

6 

4 

5 


T_TC, E_TC, Interval, SD, and Met stand for tunneling counts in temperature space, tunneling count in energy space, the 
number of MC steps between replica-exchange attempts, standard deviation, and REM based on Metropolis criterion, 
respectively. The frequency (1 & 20) of Mixed means that it was 1 MC step for random walk REM and 20 MC steps for 
designed-walk REM. 
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FIG. 2. Probability distributions of energy density at four temperatures (from left to right, 1.50, 2.10, 2.25, 2.47, and 3.15) 
obtained from (a) REM and (b) DETREM simulations including the mixed-walk simulation, and average energy density as a 
function of T from (c) REM and (d) DETREM simulations including the mixed-walk simulation and the error bars are smaller 
than the symbols. 
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